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We consider the evolutionary trajectories traced out by an infinite population undergoing 
mutation-selection dynamics in static, uncorrelated random fitness landscapes. Starting from the 
population that consists of a single genotype, the most populated genotype jumps from a local fit- 
ness maximum to another and eventually reaches the global maximum. We use a strong selection 
limit, which reduces the dynamics beyond the first time step to the competition between indepen- 
dent mutant subpopulations, to study the dynamics of this model and of a simpler one-dimensional 
model which ignores the geometry of the sequence space. We find that the fit genotypes that appear 
along a trajectory are a subset of suitably defined fitness records, and exploit several results from 
the record theory for non-identically distributed random variables. The genotypes that contribute 
to the trajectory are those records that are not bypassed by superior records arising further away 
from the initial population. Several conjectures concerning the statistics of bypassing are extracted 
from numerical simulations. In particular, for the one-dimensional model, we propose a simple re- 
lation between the bypassing probability and the dynamic exponent which describes the scaling of 
the typical evolution time with genome size. The latter can be determined exactly in terms of the 
extremal properties of the fitness distribution. 

PACS numbers: 87.10.+e, 87.23.Kg, 05.40.-a 

I. INTRODUCTION 

The episodic nature of biological evolution has provided motivation for much work on the modeling of evolutionary 
dynamics in the statistical physics community [1,2]; see [3-5] for review. Evolution displays a punctuated pattern, 
with epochs of no or slow change interspersed with bursts of (relatively) rapid activity, on various levels ranging 
from the fossil record [6,7] to experiments with microbial populations [8-11]. Punctuated behavior is also seen in 
simulations of in vitro evolution of RNA molecules [12], optimization algorithms [13,14], and artificial life [15]. It 
has been recognized for a long time that one possible scenario that is consistent with punctuated dynamics is the 
evolution of a population in a static fitness landscape with many peaks. In this picture, the two modes of evolution 
correspond to the extended periods of residence of the population at a local fitness maximum and the exploration of 
a new, higher lying peak, respectively, which entails the rapid crossing of a valley of lower fitness [16,17]. In essence, 
this is the phenomenon of quantum evolution described by the paleontologist G.G. Simpson sixty years ago [18], and 
more recently referred to in macroevolutionary theory as punctuated gradualism [6]. If the population is large, so that 
the distribution of individuals over the various phenotypes or genotypes can be modeled by a continuous field, the 
transition between fitness peaks described above displays analogies with physical processes such as quantum tunneling 
[19], variable range hopping [20] or noise-driven barrier crossing. 

Since this metastable behavior seems ubiquitous in nature, it may be worthwhile to study a simple model where this 
can be analysed in detail. A convenient mathematical framework to address this issue is the quasispecies model which 
was originally introduced to describe large populations of self-replicating macromolecules [21,22]. The quasispecies 
model is a mutation-selection model whose steady state has been studied in great detail. For various choices of fitness 
landscapes, it exhibits the phenomenon of error threshold in which beyond a critical mutation rate, the population 
delocalises over the whole sequence space. 

In this paper, we address the question of punctuated dynamics in the quasispecies model. To avoid complications 
associated with the error threshold phenomenon [4,22], we work in a strong selection limit inspired by the zero 
temperature limit of the statistical mechanics of disordered systems [23] (for a different kind of strong selection limit 
used in population genetics see e.g. [24]). In this limit, the location of the population in the space of genotypes can 
be identified with the most populated genotype at all times, and the evolutionary trajectories can be represented 
in a particularly transparent manner [25]. Since the population is located at a single genotype at any time, the 
evolutionary trajectory changes in a stepwise manner. To generate an evolutionary trajectory, a localized population 
is placed at a randomly chosen point in the space of all possible genotypes. Due to infinite population, in the next 
time step, all genotypes get occupied with nonzero population. As we describe later in detail, it turns out that the 
fit genotypes typically receive small initial population. Strong selection then reduces the problem to competition 
between these fit subpopulations struggling to overcome the poor initial conditions. 
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In the next section, we describe the quasispecies model and derive the reduced representation that arises in the 
strong selection limit. The bulk of the paper is then devoted to the analysis of the strong selection dynamics, which 
turns out to be rather closely related to the mathematical theory of records [26-29]; a relation between biological 
evolution and record statistics has been proposed previously [2,30], see also [31]. 

II. QUASISPECIES DYNAMICS IN THE STRONG SELECTION LIMIT 

The quasispecies model is defined on the space of genotypes represented by sequences a = {<ti, <tn}, where each 
of the N letters Oi is taken from an alphabet of size I > 2. The number of individuals of genotype a at time t is 
represented by a real variable Z(o~,t) that obeys the discrete time evolution equation 

Z(a, t + 1) = 5>(a' -> a) W(a') Z(a', t). (1) 

a' 

Here the fitness W{a) is defined [3] as the expected number of offspring produced by an individual carrying sequence 
a, and p(a' — > a) is the probability that genotype o is created as offspring of genotype a' due to copying errors in the 
genome. A simple choice for the latter corresponds to independent point mutations occurring with probability n per 
generation, so that p(a' -> a) = (p/(£ - l)) d ( CT »(l - ^N-d(a' ,<?) ^ where 

N 

d(a»=£(l-<^) , (2) 

i=i 

is the Hamming distance between sequences a' and a. 

A constraint of fixed population size could be enforced by dividing the right hand side of (1) by the population 
averaged fitness (W). This introduces an (inessential) nonlincarity into the problem [3], which is why we prefer to 
work with the unnormalized population variables Z(<r,t). It is clear that within the quasispecies framework, which 
requires an infinite population from the outset, the actual population size cannot play any role. 

To derive the strong selection limit [23], we write Z(a,t) = e reE ( CT >*) 5 W(a) = e KF ^ and /j, = where k is the 
inverse selective temperature [3,32] and E(a, t) and F(a) are logarithmic population and fitness variables, respectively. 
Throughout this paper we take the fitness landscape to be completely uncorrelated, which implies that the fitnesses 
F(a) are independent, identically distributed (i.i.d.) quenched random variables chosen from a distribution p(F) with 
support on interval [F m i n , F max ]. The strong selection limit corresponds to k — > oo, which yields 

E(a,t + l)=max r7 '[E(a',t) + F(a')-d(a,a')] . (3) 

Starting with an initial condition Z(a, 0) = S a CT <o) in which only a single, randomly chosen sequence has nonzero 
population, at the next time step each sequence gets seeded with the logarithmic population 

E(a,l) = F{a^)-d(a,a^). (4) 

In terms of the original variables Z(a,t), this corresponds to a population distribution that decays exponentially 
with increasing Hamming distance from . Remarkably, it turns out that for the subsequent time evolution the 
mutations are unimportant for genotypes with high fitness. Since we are primarily interested in such genotypes, the 
dynamics of the model can be approximated by allowing each genotype to reproduce with its intrinsic rate F(a) for 
t > 1. In terms of the logarithmic variables, this leads to 

E(a,t) = E(a,l) + (t-l)F(a). (5) 

Thus, (3) reduces to a problem of non-interacting sequences whose population is growing linearly in time. The 
approximation of ignoring mutations after the seeding stage was tested numerically [25] and found to be in very good 
agreement with the full model described by (3). 

A further simplification can be made by noting that due to (4), the population E(a, 1) is same for all the sequences 
that lie on a shell of constant Hamming distance d(a^°\a) away from the sequence . Since we will be interested 
in the sequence with the largest population at any instant, only the sequence with the largest fitness within a shell 
needs to be considered. Labeling a shell by its Hamming distance k = 0, 1, N from the sequence a^°\ we arrive at 
the shell model [25] in which the population E(k, t) of the fittest sequence in the shell k with a total number of ct^ 
sequences obeys 
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E(k,t) = -k + F(k)t. (6) 

Here F(k) is the maximum of at i.i.d. random variables drawn from the fitness distribution p(F); hence the F(k) arc 
independent but non- identically distributed random variables. To arrive at the simple form (6), we have redefined 
the population as E(k, t) + F(k) — F(0). The number of sequences in shell k is given by the expression 

« fe =(^)(^l) fe , (7) 

as can be seen by noting that there are (^) ways of choosing k letters at which a sequence a differs from , and 
each of these k letters can take I — 1 values. For large N, the majority of sequences is contained in a belt of width 
~ V^/V around the distance /c max = N(l — l)/£ where (7) is peaked. 

We will also consider an i.i.d. model for which the population E(k, t) evolves according to (6) but the fitnesses F(k) 
are i.i.d. random variables. This choice corresponds to a one-dimensional sequence space and we will see that our 
results are sensitive to the geometry of the sequence space. The i.i.d. model is closely related to the zero temperature 
limit of the problem of directed polymers in the presence of columnar defects studied in [33] and a version of the 
parabolic Anderson model [34]. 

Figure 1 illustrates the geometric picture of the evolutionary process (6) that emerges after these simplifications and 
approximations. At a given time t, the most populated genotype k* for which E(k*,t) = maxk{E(k,t)} leads until it 
is overtaken by a sequence k'* and so on. While the last leader is obviously located at the global fitness maximum, 
the identity of the previous leaders is non-deterministic. The inset of Fig. 1 shows the punctuated evolution of the 
leading sequence k*. Our main interest in the present work is in the statistical properties of these leadership changes 
or evolutionary jumps. Initially, the sequence with population E(0,t) leads. A shell k' can overtake the currently 
leading shell k < k! at the crossing time T(k, k') given by 

which is positive provided F(k') > F(k). In Fig. 1, only the shells k — 1,2,6,8,15 can overtake k = (and only 
k = 2,6,8, 15 can overtake k = 1, ...) since F(15) > F(8) > F(6) > F(2) > F(l) > F(0). Such a set of fitness in 
which each value is progressively larger than all the previous ones defines a sequence of records. However, in order 
to appear in the evolutionary trajectory, it is not sufficient to be a record. As shown in Fig. 1, the shell k = 2 
is bypassed by shell k = 8 since T(8, 1) = min / t > iT(fc, 1). In general, if the current leader is in shell k, then the 
next leader is in the shell k' for which the crossing time (8) is minimized [13]. Thus, in principle, the properties of 
leadership changes can be formulated in terms of the extremal statistics of the matrix of crossing times. However, as 
will be explained in more detail in Sect. VII, this minimization problem is too cumbersome to handle analytically and 
is further complicated by the presence of strong correlations among the crossing times. 

The nontrivial dynamics of the change in leadership described above is due to the competition between the initial 
population at a sequence and its fitness. As discussed above, a potential leader must be a record and hence must 
occur at a large Hamming distance away from the sequence a^' but the initial population at such a sequence is small 
due to (4). Thus, the disadvantage due to poor initial condition may be overcome by a better fitness. Further, since 
the fitness of the leader is correlated to its Hamming distance from the sequence the fitness F(k*) also shows 
punctuated evolution. 

The rest of the paper is organised as follows. In Sections III and IV, we define records and jumps precisely and 
study some of their statistical properties, both for the shell model and for the i.i.d. model. In Section V, we describe 
the dynamics of the approach to the global maximum of the fitness landscape. A conjecture concerning the probability 
of bypassing in the i.i.d. model, based on the results obtained in Sections IV and V along with numerical evidence, 
is presented in Section VI. In Section VII, we introduce and discuss a further simplification of the problem described 
by (6), which highlights the strong interdependence of the crossing times. Finally, we summarise our results and 
conclude with a brief discussion of the possible biological relevance of this work in Section VIII. 

III. RECORD STATISTICS IN SEQUENCE SPACE 

In a sequence {Xk} of random variables, an upper (or lower) record is said to occur at m if X m > Xk (or X m < Xk) 
for all k < m [26-29]. Since only the upper records are pertinent to our study, we shall henceforth refer to them as 
records. In the following subsections, we study some characteristics of these records; in particular, we find the mean 
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number of records and the typical spacing between them. The record statistics of i.i.d. variables is well studied and 
we briefly review some of the known results as they are useful for later discussion. We study the record statistics in 
the shell model, for which the random variables are not identically distributed, in some detail. An appealing general 
feature of record statistics is that the results are independent of the underlying probability distribution p(F) , which 
therefore does not need to be specified in this section. 



A. Number of records 



In this subsection, we calculate the average number 7?. s hcii of records in the shell model. For i.i.d. random variables, 
it is well known that the average number of records among A variables is 7?-iid ~ In AT for large N. To see this, 
note that the probability Piid(fc) that the fc'th random variable is a record is equal to the probability that it is the 
largest among the first k random variables. Since the location of the global maximum is uniformly distributed for 
i.i.d. variables, it follows that Pnd(k) = 1/fc for any choice of p(F). Summing Pid(fc) up to k — A yields the average 
number IZad of records to be In A. In fact, it is known that the probability distribution of the number of records is a 
Poisson distribution with mean In A [2] . 

The record statistics for non-i.i.d. variables is much less studied. A class of models in which records are obtained from 
a sequence {X k } where each X k itself is the maximum of a set of a k i-i-d. random variables was considered by Nevzorov 
[35]. Such models have been used, for instance, in an (unsuccessful) attempt to account for the frequent occurrence 
of Olympic records due to an increasing world population [36]. The i.i.d. model is a special case corresponding to 
a k = 1 for all k and the shell model is obtained by using a k given by (7). 

To analyze the Nevzorov model, we define a binary random variable Y k which takes value 1 if a record occurs in 
the fc'th set and otherwise. Since the distribution Pk{F) of the maximum of a k i.i.d. random variables is given by 
[37] 



p k (F)=a k p(F) (jf p{x)d. 




(9) 



we have 

r-F max fc-i 



Prob(F fe = 1) = / ^ dF J] qi (F) p k (F) = a \ , (10) 

Jf^ fJ- {a + ... + a k ) 

where q k (F) is the cumulative distribution corresponding to p k {F). The above equation expresses the fact that the 
fc'th record value is a global maximum amongst X)j=o a o variables and there are a k ways in which it can occur in the 
fc'th set. Further, the joint probability Prob(Y kl = l,Y k2 = 1) for k\ < k 2 is given by 

fe i-! /-F max k 2 -l 



a kl a k2 



Prob(F fel = l,Y k2 = 1) = / dF Jl qi (F) Pkl (F) / dG } \ qj(G) Pk2 (G) = j-— 

J Ftnin f£ Jf J= \f- +1 {a + ... + a kl ) {a + ... + a k2 ) 

(11) 

In a similar manner, it can be shown that Prob(Yfe 1 = 1, ■■■,Y km = 1) = Jljli ProM^fc., = 1) for any m > 2. Thus, 
the Y/c's are independent, non-identically distributed variables [27-29]. 

For the shell model, due to (7) and (10), the probability P s h e ii(k,N) that a record occurs in the fc'th shell is given 
by 

P shell (fc, N) = Prob(r fc = !)«!- (J-) -i- , a = A < Lzl = % , (12) 



1J 1 -a ' N I N 

where we have used that ( k ^ m ) / (^) ~ [a/ (1 — a)] m for k, N » 1 with k/N fixed. Since it is easier to break records 
in the beginning, the probability to find a record is near unity for k -C N. However, it vanishes beyond fc max because 
the global maximum typically occurs in the shell fc max . The average number 7^ s hcii of records can be obtained by 
simply integrating P s hcii(fc, A) over k and we find 

^ffziEiz±> w . (13) 
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Thus, 7?. s hcii increases with £ and as £ — > oo, 7£ s heii ~~ * -W because the population in each shell becomes infinite. Since 
the iVs are independent random variables with finite mean and variance, the number of records satisfies the central 
limit theorem and becomes Gaussian for large N. Specifically, for large N, the variance of the number of records is 

N N / f + 1 \ 

V she ii = ^Prob(F fe = l)[l-Prob(y fe = l)]« — — I -^-ln^-2) , (14) 

k=0 - V - / 

which is maximal for I — 5 and vanishes for large I. The ratio V s hcii/7^shcii is always considerably smaller than unity, 
which implies that the distribution is much sharper than the log-Poisson distribution obtained in the i.i.d. case. 



B. Inter- record spacings 

For the discussion of the spacings between records, it is convenient to label the records "backwards" in time, with 
r\ denoting the position of the last record (i.e., the global maximum), r 2 < r\ the penultimate record, and so on. 
In this way the pathologies associated with the fact that the expected waiting time for the next record to occur is 
infinite in the i.i.d. case [26] can be avoided. The probability P(rj) that the j'th record occurs at location rj can be 
found for the Nevzorov model in a manner analogous to (11). We obtain 

Hrj)= T TT — . (15) 

n,...,rj-i fe=0,. ..,3-1 

with N > n > ... > Tj-i > Tj and ro = N + 1. The factors on the right hand side of the above equation simply 
reflect the fact that the record at location rk+i remains the maximum until the next record occurs at > rk+i- 

For the i.i.d. model, using = 1 for all k in (15), the average location (rj) = X) r i-P( r i) can be calculated. One 
finds that the average inter-record distance Ajid(j) = (rj) — (rj+i) between the j'th and (j + l)'th record behaves as 
[28,29] 

Aii d (i)«|Qy , j = l,2,... . (16) 

A simple argument, useful for later discussion, can also be employed to obtain the above equation. The record labeled 
by j = 1) being the global maximum, is equally likely to occur anywhere between 1 and N. Thus, on average, it 
is located at N/2. Similarly, the record labeled by j = 2 is a global maximum in the range [l,ri) with uniformly 
distributed location, which gives the average location (r 2 ) = (1/2) (n) = N/4. Repeating this argument, we obtain 
the result in (16). 

For the shell model, since the most likely position of the global maximum is in the shell with the largest number 
ctk of sequences, we have (r\) = k max . For sake of simplicity, we consider binary sequences (I = 2) in the following 
but the scaling behavior obtained below holds for I > 2 as well. We find that for j > 2, the average location (rj) of 
the j'th record is given by 



<rj> = <ri> - 2 {v^J J-J^- 1 L, dx >- 2 crf^y - L dxi rf&D ' j - 2 ' (17) 

where we have used Eqs.(Al)-(A3) and performed a Gaussian integral. The average location (r 2 ) of the second record 
is given by (r 2 ) = (ri) — 2.0064V1V /ir^/2. Thus, the second record (and in fact, j'th record for j of order unity) can 
be found within 0(y/N) distance of the global maximum since at has width <~ \fN about fc max . For j > 2, after 
repeated integration by parts, (17) can be rewritten as 

, < , < 1 fN (ln2V'- 2 - m 

fa> = <n> - - V T E G(m) > * > 2 ' (18) 

m=l VJ ' 



where 



G(m) = (-1)™ / d X ^L. (Inerfc(x)r-(ln2r (1Q) 
, erfc(x) ml 
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As outlined in Appendix A, the integral G(m) can be estimated by the saddle point method and we find G(m) w 
for large m. Since the leading contribution to the sum in (18) comes from the m = j — 2 term, we have 



~ / 1 V 

A she ii(j) « J — , J > 1 • (20) 

Thus, while the inter-record spacing decays exponentially with j for the i.i.d. model, it falls as a power law for the 
shell model. The spacing between the first few records [j = 0(1)] is of order >/~N, while for the bulk of the records 
with j = 0(A) the spacing is of order unity; this is consistent with the vanishing of the record occurrence probability 
(12) near k = fc max . 



IV. JUMP STATISTICS 



As we described in Sect. II, in our model, evolutionary jumps are a subset of records and if a jump occurs at k', 
the next jump is said to occur at k > k' if (i) F{k) is a record (ii) the overtaking time T(k,k') = mm,->fc/ {T(j, k')}. 
By convention, the first jump and the first record occurs at k = 0. Due to the second condition, some of the records 
can get bypassed and fail to appear in the set of jumps. In this section, we find the mean number of jumps and the 
inter-jump spacing for both the i.i.d and the shell model. In contrast to the properties of records, the statistics of 
jumps depends explicitly on the fitness distribution p(F) and we will consider distributions for which the tail behavior 
corresponds to the three universality classes of standard extreme value theory [38] and which also appear in a very 
similar form in the theory of records [27-29] . 



A. Mean number of jumps 

We begin by discussing numerical results for the average number of jumps in the i.i.d. model. It was found in [25] 
that the average number Jlid of jumps grows as /31n A where the prefactor /3 < 1, and was conjectured to be 

f 1/2 , p(F) ~ e~ F 

/?* { (<J-l)/(25-l) , p(F) ~ F' 1 - 5 , 6>1 . (21) 

[ (2 + i/)/(3 + 2u) , p(F) ~ (F max - FY, v>-l 

Figure 2 shows J\\<± increasing linearly with In A and slope (3 for some distributions in accordance with (21). Thus, 
in the i.i.d. model the jumps can be viewed as "diluted" records, in the sense that the mean number of records and 
jumps differ only up to a prefactor and the probability Piid(fc, N) that a jump occurs at k is given by p/k. In this 
picture, the probability for a given record to be bypassed is simply 1 — 0. However, bypassing is not completely 
random, as the variance of the number of jumps is found consistently to be smaller than the mean. This implies 
a certain amount of "anti-bunching" among the jumps, which can also be detected by the direct measurement of 
correlation C;id(fci, ^2) = Pud(ki, £2) — Pud(ki)Piid(k2) where Pnd(ki, k?) is the joint distribution of having a jump at 
k\ and &2, as shown in the inset of Fig. 4. Some further discussion of the conjecture (21) will be provided in Sect. VI. 

Estimates for the mean number of jumps in the shell model were also given in [25] , but the range of sequence lengths 
was too limited to allow for a definite statement. Here we present the results of our simulations for large values of A 
obtained using the approximation described below. Since the shell fitness F{k) is the maximum of i.i.d. random 
variables drawn from the distribution p(F), it can be obtained from a uniform random variable u using the relation 

r F(k) 

J p(x)dx = u 1/a ». (22) 

However, the binomial coefficient (^) increases exponentially with A for large k and the distribution of u 1 / Qfc ap- 
proaches a delta- function centred at unity for k 3> 1 making it difficult to determine the distribution pk{F) of F(k) 
accurately when N is large. For the exponential fitness distribution, the relation (22) becomes 

F(k) = - ln(l - e lnu/ak ) ~ - Inpn(l/«)] + In(a fc ). (23) 

Since the last expression only involves the logarithms of binomial coefficients, F(k) can be easily generated up to 
large values of N. While a similar approximation can be employed for other distributions with unbounded tails, we 
have not been able to obtain reliable results for bounded distributions. 
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In Fig. 3, we show simulation results for the probability P s he\i(k, N) that a jump occurs in the fc'th shell, for a 
binary alphabet (£ = 2) and two different values of N. The data collapses onto the scaling form 

P shen (k, N) « N-Wfth/N) , x < = ^ , (24) 

where the scaling function behaves as f(x) ~ x -1 / 2 for small x. This has the interesting consequence that P s h e u ~ 
l/y/k is independent of TV, and hence the number of jumps grows as \fk for k <C fc max . The scaling function appears 
to be independent of the alphabet size £, which therefore only changes the cutoff at k = fc max [31]. Thus, the average 
number J7" s heii of jumps obtained by integrating P s heii(&> N) over k grows with N as 

JshcU-yJ^pE . (25) 

Our simulations indicate that the above dependence of i7 s heii on N is also true for normal-distributed fitness and we 
expect it to hold for all distributions decaying more rapidly than any power law [31]. Further, the distribution of the 
number of jumps is a Gaussian (as in the case of records) with both mean and variance scaling as y/~N . 

For the power law case, we find that J^eii — » 1 for large N. As we shall see in Section V, the globally fittest 
sequence takes over the leadership in a time of order unity in this case, which explains the above behavior of J s h e u- 
In summary, unlike the i.i.d. model, most of the records are bypassed in the shell model both for exponential and 
power law distributions. 



B. Inter-jump spacings 



We now turn to a discussion of the inter-jump spacings A(j) defined in analogy to the inter-record spacings. We 
denote by Sj the position of the j'th jump, with j — 1 referring to the last jump, which is also the last record (the 
global maximum), and define the inter-jump spacing as A(j) = (sj) — (sj+i) where (sj) is the average location of the 
j'th jump. For the i.i.d. model, an approximate calculation of Aiid(j) can be carried out by assuming jumps to be 
randomly diluted records. The position Sj of the j'th jump equals the position of the fc'th record, where k > j 
because of the possibility of bypassing. If the record k + 1 is not bypassed, then sj + i = r k+ i and Sj + i — (l/2)sj on 
average due to the argument given after (16) for the inter-record spacings. In the diluted record picture, this is true 
with probability (3. With probability (3(1 — (3), record fc + 1 is bypassed and s J+ i = so that Sj+i = (l/4)sj on 
the average. Similarly, with probability (3(1 — (3) 1 , I records are bypassed and the ratio Sj + i/sj = (l/2)' +1 on average. 
Summing over all possibilities, we find that (sj+i) = b(sj) with 

oo „ 

6 = £ 2 -(m)^ (1 _^ = _A_. (26) 

Taking into account that the sum over all inter-jump spacings adds up to (s\) = (ri) = N/2, we obtain 

AiidO')«f (I"*)* 7 '" 1 - i = 1,2,... . (27) 

The numerical data shown in Fig. 4 supports the general form of this expression, but with the coefficient b ^ (3/2 
rather than (3/(1 + (3). This is another indication of correlations that are not accounted for in the random dilution 
picture. 

For the shell model, as shown in Fig. 3, the data for A s h e ii(i) for various N collapses onto a monotonically decreasing 
curve if we assume A s h e ii(j) to be of the scaling form 

A shen (j)~VNh(j/VN) . (28) 

The scaling with v^/V follows naturally from the scaling form (24) for the jump occurrence probability. The latter 
implies that the density of jumps on the fc-axis is of order l/y/N for finite k/N, hence the spacing is of order y/N, 
and the argument of the scaling function is j/ \/N because the total number of jumps is also of order \/N. However, 
the tail of the scaling function h does not obey the scaling form (28) and approaches zero with increasing N. Thus, 
while for finite k/N, the jumps are roughly equally spaced with spacing \/N, the spacing is of 0(N S ) with s < 1/2 
forfc/TV^O. 
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V. APPROACH TO THE GLOBAL FITNESS MAXIMUM 



So far, we have discussed the statistics on the k* axis of the inset of Fig. 1 and now we focus on the temporal 
statistics. As explained in Section II, a fit sequence k* leads till it is overtaken by an even fitter one and eventually 
the globally fittest sequence emerges as the leader at typical time T*. In this section, we find this typical time T* 
required for the population to reach the global maximum of the fitness landscape and the distribution of the evolution 
times to jump from one local maximum to another. 



A. Dynamic scaling 

The location k*(t) of the most populated sequence at time t for which the logarithmic population is E(k*,t) = 
maxfc{_E(fc, t)} increases with time till the global maximum is reached. In Appendix B, the distribution P t {k*) is 
explicitly calculated for the i.i.d. model in the limit of infinite genome size (N — > oo). This distribution is usually of 
the scaling form 

P t (k*) =t- 1 ' z <S>{k* /t 1 '*) , (29) 
where the scaling function $ depends on the underlying fitness distribution and the dynamic exponent z is given by 

1 , p(F)~e- F 

(2 + + v) , p(F) ~ 0F max - F) v , v>-\ , (30) 
(5-l)/8 , p(F)^F- 1 - s , 6>1 



for the three classes of fitness distributions introduced in Sect. IV. The corresponding behavior k*(t) ~ t x / z of the 
mean location has been derived previously using Flory-type arguments [20,33] and can also be seen using (29). As 
we have already discussed, the global maximum is reached when k*(t) w N/2, which defines a total evolution time 
T* <~ N z . One thus expects a scaling form 

k*(t,N) ss t 1/z p(t/T*) , (31) 

where the scaling function <p(x) is a constant for x -C 1 and decays as x~ x l z for x ^> 1 [39]. 

An alternative approach [23,25] to estimating the typical time T* required by the population to reach the global 
maximum starts from the observation that T* is of the order of the time T\ at which the globally fittest sequence at 
typical location r\ = s\ and fitness F{r\) overtakes the penultimate leader with respective quantities s 2 and F(s2) 
(refer Sect.IVB). From (8), we have 

Tl = F(s7) = F(s 2 ) ' (32) 

where the inter-jump spacing in the numerator is given by (27). The estimation of the denominator involves a subtlety 
in previous works [23,25], it was assumed that F(si) — F(s 2 ) is of the order of the fitness gap e, which was defined as 
the difference between the values of the global maximum and the second largest fitness of the fitness landscape. Because 
the second largest fitness does not necessarily appear in the record sequence, e is only a lower bound on F(ri) — F(r 2 ), 
which in turn is clearly a lower bound on the fitness difference of interest, e < F(n) — F(r2) < F(si) — F(s 2 ). However, 
our explicit calculations show that F(ri) — F(r 2 ) is of the same order as e; moreover, at least for the i.i.d. model, we 
know that at most a few records are bypassed between si and s 2 , and hence the assumption that F(s\) — F(s 2 ) <~ e 
seems justified. 

The calculation of the distribution of fitness gap e is a standard exercise in extreme value statistics [37]. In a 
system with a total number S of sequences, the typical value of the fitness gap increases as S 1 ^ 6 for the unbounded 
power law distribution, decreases as S^ 1 ^ 1+ ' J " > for the bounded distribution, and is of order unity for the exponential 
distribution. For the i.i.d. model, using si — s 2 ~ N [see (27)] and S — N in (32), we recover T\ ~ T* <~ N z with z 
given in (30). The other fitness distributions can be treated in a similar manner. 

For the shell model, S — £ N and due to (28), the numerator is of order \/N so that T* <~ \/N for the case of 
exponentially distributed fitness. Presumably, the time Tj at which the j'th jump occurs is also of order y/N for 
j ~ 0(1)- Since the total number of jumps is of the same order, it follows that initially there are many, quick 
jumps followed by few jumps that take 0(VN) time. This result agrees qualitatively with that seen in experiments 
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(discussed later) concerning the pace of evolution which is initially rapid and later slows down considerably. For the 
bounded distributions, T* increases exponentially with N whereas it decreases exponentially for the fat-tailed power 
law distributions. The latter result implies that, for large N, the global maximum takes over in a single time step, 
which explains why the mean number of jumps tends to unity for the power law distributions (see Sect. IV A). 

B. Universal tails of the evolution time distribution 

In the last subsection, we found the typical time T* to reach the global maximum and now we consider the 
distribution P(T\,N) of the time T\ at which the final jump occurs. For the i.i.d. model, since the typical T\ also 
grows as N z , we may expect the normalised distribution P(Ti,N) to be of the scaling form P(T\,N) N~ z gi(Ti/N z ). 
In general, for j ~ 0(1), the distribution P(Tj,N) of the time Tj at which the j'th jump occurs is of the scaling form 

P(T j ,N)^N- z g j (T j /N z ) . (33) 

Although the dynamic exponent z depends on the underlying fitness distribution, we shall now show that the tail of 
the distribution P(Tj,N) is universal. The events contributing to large T\ are the ones for which F(s\) — F(s 2 ) is 
small; in these cases we expect the general bound F(s\) — F(s 2 ) > F(r\) — F(r 2 ) to be saturated, i.e. the second 
record is not bypassed and s 2 = r 2 . Thus, using (32), we obtain 



P(T U N) 



dTi 



Prob( £l = A iid (l)/T!) « Prob( £l = 0) , (34) 



for large T 1; where €\ — F{r\) — F(r 2 ). The probability distribution of €\ can be obtained along the lines of the 
derivation of the distribution of the fitness gap e in [23], and it is found that the probability for a near- vanishing 
difference between two successive record values is nonzero for any p(F). We conclude that P(T\,N) has a power law 
tail with exponent —2 for any underlying fitness distribution [25]. This is an example of the generation of a power 
law through a change of variables (from ei to T\ ~ 1/ei) as described in [38]. 

Similarly, the events contributing to the tail of P(T 2 , N) are the ones in which the record at r 2 is the penultimate 
leader and the record at r% is the leader previous to it. Thus, we demand that none of these two records should be 
bypassed; in particular, r 2 should not be bypassed, which requires that T(n,r 2 ) > T(r 2 ,r 3 ). This condition can be 
written as t\ < Ce 2 , where e 2 = F(r 2 ) — F(r 3 ) and C is the average value of (n — r 2 )/(r 2 — r 3 ), a number of order 
unity. Thus we obtain 

A- ,(2) r CAnd{2)/T2 CA-A2) 2 
P(T 2 , N) « ^£2 / d€l Prob( ei , e 2) N) « ^ ! Prob( £l - 0, e 2 = 0, N). (35) 

J 2 JO i 1 

Since Prob(ei = 0, e 2 = 0, TV) can be shown to be nonzero, we conclude that P(T 2 ,N) - T 2 ~ 3 for large T 2 . Extending 
the above arguments in a similar fashion to the next evolution times, we find that the scaling functions in (33) behave 
as gj(x) <~ x~ x ~i for x 3> 1. Interestingly, this implies that the expected time (Tj) is finite for j > 2 and infinite 
for j = 1. In Fig. 5, the prediction P(Tj,N) <~ T~ X ~' J for j = 1,2,3 is compared with data obtained using Monte 
Carlo simulations for the i.i.d. model. The behavior of the universal tails of P(Tj) discussed above is true for the 
shell model as well. 

VI. BYPASSING PROBABILITY AND DYNAMIC EXPONENT: A CONJECTURE 

As we have seen in the previous sections, the jump statistics are not analytically tractable due to the constraint of 
minimal overtaking time. For the i.i.d. model, the jumps differ from the records only up to a prefactor (3 conjectured 
to be given by (21). Comparing the expressions in Eqs.(21) and (30), we observe that the bypassing probability 1 — (3 
appears to be related to the dynamic exponent z by the following universal relation 

0=(l-0)z . (36) 

A derivation of this relation (which eludes us so far) would constitute a proof of the conjecture (21). 

Interestingly, the relation (36) can be interpreted in terms of a kind of duality between the k- and the i-axis of the 
inset of Fig. 1. So far we have identified each jump with the position on the .E-axis where the line that takes over the 
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leadership when the jump occurs originates; but we may just as well identify the jump with the corresponding crossing 
time T(k, k') at which the leadership shifts from k' to k. Clearly there is a one-to-one correspondence between the 
jumps defined on two axes and the average number J( iA of jumps on the time axis is equal to Jud discussed in earlier 
sections. Thus, J( {d — ~ /31n N and since the typical time T* to reach the global maximum scales as N z , we have 
J( id « 0' InT* with 0' = 1 — due to the conjecture (36). This leads us to expect that the probability Pad(t, N) that 
a jump occurs at time t decays as 0'/t for t <C N z and as 1/t 2 for t » N z ; the latter behavior is the universal t~ 2 
tail explained in Section V. We conclude that the sum of the jump probabilities along the k- and the i-axes should 
sum up to a universal function, i.e. 

P iid (t = X, N) + P iid (k = X,N) = l/X , (37) 

for any choice of fitness distribution. The numerical evidence supporting this claim is shown in Fig. 6. Furthermore, in 
analogy to the jump spacing A;id(j) along the fc-axis, one can also consider the quantity A' iid (j) = (Tj) — (Tj + \) which 
is the spacing between the successive jumps on the time-axis. Replacing N by T* and b = 0/2 by 0' /2 = (1 — (3)/2 
in (27), we expect 

A(id0-)~^(^) 3 1 . i = l,2,... . (38) 

Numerical results consistent with this expression are shown in the inset of Fig. 6 for some distributions. The deviations 
seen in the data for the first jump (j = 1) reflect the fact that, because of the 1/Tf-tail derived in Eq.(34), the average 
of T\ is not defined and hence grows with the number of disorder realizations. 



VII. A MODEL BASED ON RECORD TIMES 



In this section, we introduce a further simplification of the i.i.d. model. As we have discussed already, a sequence k 
may occur in the set of jumps provided F(k) is a record. Thus, it is sufficient to consider only the subset of sequences 
whose fitness is a record. Here it is convenient to label the records forward in time, so we denote by Rj the location 
of the j'th record with j = 1 labeling the first record (i?i = 1 by convention), i?2 > Ri the second record, and so 
on. Note that there are two sources of randomness in the problem - one arising from the record locations Rj and 
the other due to record values F(Rj). For exponentially distributed fitness, it is known that the differences between 
successive record values are independent and exponentially distributed random variables [40]. Thus, the fitness of two 
successive records differs by unity on average. These considerations allow us to eliminate the randomness associated 
with the record values by replacing the i.i.d. model of (6) with exponentially distributed fitness by a simpler model 
for which the population evolves as 

E(j,t) = -Rj+jt. (39) 

Like in the original i.i.d. model, we find numerically that the average number J of jumps grows logarithmically with 
N with a prefactor (3 w 0.63 (see Fig. 7). This is distinctly different from the value (3 w 1/2 found for the i.i.d. model 
with exponential fitness distribution, indicating that the randomness in the record values is relevant. Somewhat 
surprisingly, in contrast to the conjectured values of (3 for the full i.i.d. problem given in (21), (3 does not seem to be 
a simple rational number. 

Our primary motivation for introducing this simplified model is to gain further insight into the mathematical 
structure of bypassing. For the model defined by (39), the crossing time T(j,j') at which the line associated with the 
j'th record overtakes that associated with record j' is given by 

/•:./•/: ^-4- (40) 
3-3 

Then the probability $2 that the second record is not bypassed can be written as 

i?3 — 1 R4 — 1 



fa = Prob 



R2 — 1 = min i? 2 — 1, 



(41) 



The evaluation of the condition on the right hand using the joint probability distribution for the record times [27-29] 
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Prob( fl2 , ._, H.) . (ft _ 1)(R3 _ im _ 1} ^ _ ^ ,42) 

is clearly a difficult task. An upper bound on /3 2 is obtained by requiring that the record at R 2 is not bypassed by 
the one at i? 3 , i.e. that T (2, 1) < T(3, 1). This gives 

oo 

& < Prob(i? 3 > 2R 2 - 1) = V - — = 2 (1 - In 2) w 0.613706 . (43) 

£"^ 2 Jt2 - 1 Zti,2 - 1 

In our simulations on a large system, we find /3 2 ~ 0.600786, showing that bypassing of i? 2 by the records beyond R 3 
is rather unimportant. 

Consider next the behavior of the crossing times (40) when j and f are large. Williams [41] has shown that the 
sequence of record times can be generated from the recursion relation [26-29] 

R j+1 = [e x 'R j ] + l, (44) 

where the Xj are independent, exponentially distributed random variables with mean one, and [a] is the integer part 
of a. For large j, the integer constraint can be ignored, and hence 



r 



exp ( £ A, ] - 1 

i=j'+i 



(45) 



Recalling that the choice of the next non-bypassed record involves finding the minimum among all crossing times 
T(j,f) with j > j', we see that that the current location Rj> cancels in the comparison between two such crossing 
times. The problem thus acquires translation invariance in the record space, in the sense that the position of the next 
non-bypassed record depends only on j — j' and on the random variables X; t associated with the records between j' 
and j. It is therefore plausible that the bypassing probability tends to a constant for large j, and one is tempted to 
describe the process by a Markov chain on the set of records with the transition probability 

P rj = Prob[f (j,f) = mint (n,j% (46) 

This is the conditional probability that the next jump occurs at j, given that the preceding jump was at j', averaged 
over all realizations of record times. Using the representation (45), the Pjji are manifestly translationally invariant 
for large depending only on j — j' . Even in the asymptotic limit in which the expression (45) can be used, the 
evaluation of (46) is cumbersome, but an analytic upper bound on Pj,j+i can be obtained along the lines of (43). In 
the limit of large j, we can write 

Pjj+i < Prob[T (j + 1, j) < f(j + 2, j)] = Prob[e^+ 1 + e~ x ' > 2] = In 2, (47) 

where Xj and Xj+i are the independent exponential random variables used in the representation (44). The numerical 
cvalution of (46) yields Pj,j+i ~ 0.669, Pj,j+2 ~ 0.225 and Pj,j+3 ~ 0.075, indicating a roughly exponential decay 
of the transition probability. From the Pj>j, the mean density of jumps (non-bypassed records) can be computed 
according to 

^Markov = ^XjnP, ij+ „^ w 0.676, (48) 

which is significantly larger than the direct numerical estimate (3 « 0.63 (see Fig. 7). This shows that the transition 
probability (46) is not an exact representation of the process. The reason is that (46) is an annealed average, whereas 
in the full problem the record times (or, equivalently, the exponential random variables in (44)) must be treated as 
quenched: Minimizing the crossing times T(j,f) for a given j' involves, in principle, all j > j', and the same set of 
random variables is used every time this minimization is repeated for different j' . 

We thus have to conclude that the range of attainable analytic results, even for the simplified problem (39), is very 
limited. An extension of the bound (43) to the full i.i.d. model should be feasible using the representation of the joint 
distribution of record times and record values through a Markov chain [28,29]; however, as such a bound is unlikely 
to provide much insight into the conjectured relation (36), we have not pursued this approach. 
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VIII. CONCLUSIONS 



In this article, we characterised the evolutionary trajectories traced out by a quasispecies population in an un- 
corrected rugged fitness landscape. These trajectories approach the global fitness maximum through a sequence of 
jumps which mark a change in the identity of the most populated genotype. The statistics of these evolutionary 
jumps was studied mainly numerically. However, useful insights were provided by a study of record statistics which 
could be handled analytically. It was found that the jump statistics are qualitatively similar to records, but there are 
quantitative differences because, as shown in Fig. 1, a record breaking genotype can be bypassed by a superior one 
before it can acquire dominance in the population (i.e. qualify to be a jump). 

The statistics of records and jumps depends strongly on the geometry of the space of genotypes. The natural setting 
for genotype evolution is the Hamming space of sequences of fixed length N. However, computational effort could be 
greatly reduced by lumping together the sequences within a shell of constant Hamming distance with respect to the 
initial population [25]. Complementary to this shell model, we also considered a model of i.i.d. shell fitnesses, which 
corresponds effectively to a one-dimensional sequence space. While for the i.i.d. model, the average number of jumps 
differs from the number of records only through the prcfactor (3 of the logarithm of N, for the shell model the ratio 
of the two numbers J7sheii/7^shcii — ► as N — > oo. For fat-tailed fitness distributions the evolutionary trajectories in 
the shell model may even degenerate, in the sense that the global fitness maximum is reached in a single step. For 
distributions decaying faster than a power law, like the exponential and normal distributions, we find numerically 
that J7 s hcii ~ V^V ; an analytic understanding of this result would be very desirable. 

A universal feature of the evolutionary trajectories, which is independent of the geometry of genotype space as well 
as of the fitness distribution, is the hierarchy of power law tails for the distributions of the times at which the jumps 
occur. In particular, the T~ 2 -tail for the total evolution time implies that the average of T is infinite. The dependence 
of the typical evolution time on the size of the sequence space can be characterized by a dynamic exponent z, which 
was obtained exactly in terms of the extremal properties of the fitness distribution. On a mathematical level, perhaps 
the most intriguing result of this work is the conjectured relation (36) between (3 and z for the i.i.d. model. It would 
be extremely interesting to understand how such a simple relation arises from the properties of the matrix of crossing 
times (8); however, in view of the difficulties encountered even in the analysis of the simplified problem (39), we do 
not see at present how further progress in this direction can be achieved. 

In this article, we worked in the limit of infinite population and strong selection. However, we expect our results 
to hold, at least qualitatively, for finite selective temperature as well. The jumps appear instantaneous in the strong 
selection limit; at finite selective temperature they take a finite amount of time in which the peak associated with the 
new leader catches up with the currently dominating peak and the population distribution briefly becomes bimodal 
[13,19,33]. Although the quasispecies theory, which works in the infinite population limit, is inadequate to address 
the fluctuation effects that become important when small mutant populations cross a fitness valley [42], for the 
related problem of episodic behavior in evolutionary computation [14], some of the finite population behavior has 
been understood on the basis of the infinite population limit. Thus we expect our investigation to give some insight 
into the behavior of finite populations in rugged fitness landscapes [43,44] as well. 

We close with some remarks on the applicability of the present work to the evolution of biological populations. A 
realization of asexual mutation-selection dynamics in a static fitness landscape that is believed to be quite rugged 
is provided by the long-term experiments on populations of Escherichia coli carried out by Lenski and collaborators 
[8,9,11]. These experiments show evidence for punctuated behavior both in the fitness and in the morphological 
features (such as cell size) of the evolving populations, which is attributed to the emergence and fixation of beneficial 
mutations. Fixation implies that a mutation which is initially present only in a single individual is inherited by a 
growing fraction of the population and eventually acquires dominance. 

At least on a qualitative level, this process corresponds in our model to that in which a subpopulation residing 
in a distant shell of sequence space takes over the leadership of the population, as described by (6) and illustrated 
in Fig. 1. The bypassing of one subpopulation by another is analogous to the phenomenon of clonal interference, in 
which one beneficial mutation is superseded by another one before reaching fixation [45] . The key difference between 
our model and the behavior of real asexual evolution is that in the latter case beneficial mutants arise through the 
stochastic search of a finite population in an immensely large sequence space, while in the quasispecies model all 
possible mutants are present (in extremely small numbers) after the first time step. It is, therefore, mandatory to 
include the stochastic finite population dynamics into the model. Nevertheless, some features of the competition 
between the mutant populations (however they may have arisen) could well survive in a more complete, realistic 
treatment. 
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APPENDIX A: INTER-RECORD SPACING FOR THE SHELL MODEL 



Using the Stirling's formula r! w \/2nr(r/e) r in the binomial coefficient (^) for r, N 3> 1 with r/7V fixed and 
expanding (^) about its maximum at r = N/2, we have 



1 

2~N 



This expression can be used to show that 







(AT/2 - r) 2 ~ 


0- 




7V/2 



1 y 



2 iV 



r=0 



-erfc(a) 




erfc(a) 



7rN 



exp(— a 2 ) 



(Al) 

(A2) 
(A3) 



by approximating the sum on the left hand side by an integral. In the above expressions, a = (N/2 — y)/y/N/2 
and erfc(x) = (2/y/n) J°° dt e~* is the complementary error function. One can derive (17) for the average location 
( r j) = '52 r jP{ r j) where P(rj) is given by (15) using Eqs.(Al)-(A3) and replacing the sums by integrals. 
The integral G(m) in Section IIIB can be estimated by saddle point method as follows. We have 



G(m) = (-1)' 



e~ 2x 

dx 

-oo erfc(a;) 



(lnerfc(x)) m - (ln2) r 



i-iy 



dx 



-Q(x) 



(A4) 



where Q(x) w x 2 — In x (2m + 1 



2 ) for large m. Here we have used that crfc(x) w e x /(y/irx) for x 3> 1. 



Expanding Q(x) about its minimum x ~ m — ln-^/m up to 0(.t — a; ) 2 and doing the Gaussian integral, we obtain 
G(m) s=s yJnm/2 for large m. 



APPENDIX B: DISTRIBUTION OF THE MOST POPULATED SEQUENCE 

Our goal is to compute the probability Pt(k*~) that the sequence k* has the maximum population at time t in the 
i.i.d. model. This distribution is given by 

Pt(k*)= r max(fe dE pf\E) \{ q [ k) (E) , (Bl) 

where 

p { t k) (E)^t- 1 p[(E + k)/t} (B2) 

is the distribution of E(k, t) obtained from the fitness distribution p(F) via the variable transformation (6). The limits 
of the support of p [k) (E) are E min (k,t) = -k + F min t and £ max (M) = -k + F max t, and q { t k) (E) = J^dE p {k) (E) 
is the corresponding cumulative distribution for E > E m - ln and zero otherwise. In the following, wc show that the 
distribution Pt(k*) is of the scaling form Pt(k*) rts t~ x / z $(fc* /t 1 ^) for various choices of fitness distribution, 
(i) p(F) = e- F : 

p t (k*) = - dx e -(*+ fe *)/* rr (i - e -(-+fe)/*) „ i e - fe */* , (B3) 
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where the last expression is obtained by exponentiating the product and evaluating the resulting sum as an integral 
for an infinite system. Thus in this case z — 1 and the scaling function is given by &i(y) — e~ v . 
(ii) p{F) = \l + u){l-FY,u>-V. 



Pt{k*) 



l + v 



l + v 
t 

1 



t-k* 



dx 1 



r-t-k" 

Jo 



dx 1 



k* 



x + k* 
x + k* 



n 

k=£k* 

exp 



1 - I I 

t 



2 + 



;0-f) 



2+1/ 



where z = (2 + u)/(l + v) and the scaling function $2(2/) is given by 

1 + f 



$2(2/) 



(2 + !/)!/* J y 2+, /{2+v 



dx 



1 / z (y+{{2 + V )xY'^) 1 ' . 



(iii) =&F~ 1 -\ 5>l: 



x + k* 



1+5 



c />00 

p t (fc*) = ^ y ^ 

where z = (5 — l)/5 and the scaling function $3(2/) is given by 



n 



x + k 



Jo 



dx 



e -x x l/(S-l) 



(1 + ((S -l)x) 1 /( s - 1 )y) 1 + s ' 



(B4) 

(B5) 
(B6) 

(B7) 



(B8) 
(B9) 
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FIG. 1. Time evolution of the logarithmic population variable E(k,t) — — k + F(k)t. The fitnesses F(k) of the sequences 
corresponding to thin green lines are not records whereas those corresponding to bold red lines are. The lines appearing in the 
upper envelope define the most populated sequence k*; an evolutionary jump occurs when two such lines cross. Note that the 
sequences k — 2 and k = 6, although constituting fitness records, are bypassed as they do not satisfy the minimum crossing 
time constraint. The inset shows the punctuated evolution of the most populated sequence k* as a function of time. 
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FIG. 2. Mean number J7iid of jumps for the i.i.d. model for various fitness distributions. The lines are the best fits to the 
functional form J7iid = (UlnN + constant. 
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FIG. 4. Semi-log plot for inter-jump spacing Am(J) for the i.i.d. model in accordance with (27) with slope b = (3/2. Inset: 
Correlation Cijd(fei, k 2 ) vs. k 2 — ki for two fixed values of ki to show correlations between jumps in the i.i.d. model. 
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FIG. 5. Log-log plot for the ratios P{T 2 )/P(T 1 ) and P(T 3 )/P(T 2 ) of evolution time distributions for the i.i.d. model with 
p(F) — e~ F . The ratios decay as 1/T for large T, as the line with slope —1 indicates. 
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FIG. 6. Distribution Pi^{X, N) — Pad(t = X,N) + Pad(k = X,N) vs. X for various distributions, illustrating the duality 
between the jump processes along the k- and t-axes. The solid line has slope —1. Inset: Semi-log plot for the inter-jump spacing 
Aj id (j) along the time axis as a function of j for the i.i.d. model. The slope (1 — fl)/2 is consistent with the conjecture (38). 
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